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Chapter 1 

Geometric Integration 



In this chapter, we will discuss the basic issues about Geometric Integration, giving a 
concrete motivation to look for energy- conserving methods, for the efficient numerical so- 
lution of Hamiltonian problems. In particular we will focus on the basic idea Hamiltonian 
Boundary Value Methods (HBVMs) relies on, i.e., the definition of discrete line integrals. 
The material of tis chapter is based on references [HI SH HOI IHl IZ] • 



1.1 Introduction 



The numerical solution of conservative problems is an active field of investigation dealing 
with the geometrical properties of the discrete vector field induced by numerical methods. 
The final goal is to reproduce, in the discrete setting, a number of geometrical properties 
shared by the original continuous problem. Because of this reason, it has become custom- 
ary to refer to this field of investigation as geometric integration, even though this concept 
can be led back to the early work of G. Dahlquist on differential equations, aimed at re- 
producing the asymptotic stability of equilibria for the trajectories defined by a numerical 
method, according to the well-known linear stability analysis (see, e.g., |25j). 

In particular, we shall deal with the numerical solution of Hamiltonian problems, which 
are encountered in many real-life applications, ranging from the nano-scale of molecular 
dynamics, to the macro-scale of celestial mechanics. Such problems have the following 
general form, 

y' = JVH{y), y{0) = y, E M^™, (1.1) 

where = — J = is a constant, orthogonal and skew-symmetric matrix, usually 
given by 

j=( i \ (1.2) 



-/ 

(here I is the identity matrix of dimension m). In such a case, we speak about a problem 
in canonical form. The scalar function H{y) is the Hamiltonian of the problem and its 
value is constant during the motion, namely 

H{y{t))^H{y,), Vt > 0, 



for the solution of (1.1). Indeed, one has: 

^H{y{t)) = VH{y{t)fy\t) = VH{y{t)fJVH{y{t)) = 0, Vt > 0. (1.3) 

Often, the Hamiltonian H is also called the energy, since for isolated mechanical systems 
it has the physical meaning of total energy. Consequently, energy conservation is an 
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important feature in the simulation of such problems. The state vector of a Hamiltonian 
system splits in two m-length components 

^( q 

y 

where q and p are the vectors of generalized positions and momenta, respectively. Conse- 
quently, (1.1)-(1.2) becomes 

q = VpH{q, p), p = -VqH{q, p) . 

Depending on the case, we shall use both notations. 

Another important feature of Hamiltonian dynamical systems is that they posses a 
symplectic structure. To introduce this property we need a copule of ingredients: 

- The flow of the system: it is the map acting on the phase space M^"* as 



where y{t) is the solution at time t of (1.1) originating from the initial condition y^. 
Differentiating both sides of (1.1) by ?/o and observing that 

dy{t) d(f)t{yo) . 
dyo dyo 

we see that the Jacobian matrix of the flow 0^ is the solution of the variational 
equation associated with (1.1), namely 

= JV'H{yit))A{t), A{0) = I, (1.4) 

where V'^H{y) is the Hessian matrix of H{y). 

- The definition of a symplectic transformation: a map u = {q,p) G M^™ u{q,p)M.'^"^ 
is said symplectic if its Jacobian matrix u'{q,p) G ]R2mx2m jg ^ symplectic matrix, 
that is 

u'{q,p)'^Ju{q,p) = J, for all q,p e M™. 

That said, it is not difficult to prove that, under regularity assumptions on H{q,p), the 
flow associated to a Hamiltonian system is symplectic. Indeed, setting 

dyo 

and considering ( |1.4[ ), on has that 

= {A{tfV^H{y{t))J^JA{t)) + {A{tfjJV^H{y{t))A{t)) = 0. 

Therefore 

A{tfjA{t) = A{OfjA{0) = J. 

The converse of the above property is also true: if the flow associated with a dynamical 
system y = f{y) defined on M^™ is symplectic then necessarily f{y) = J'WH{y) for a 
suitable scalar function H{y). Consequently, conservation of H{y) follows, by virtue of 
(0. 

Symplecticity has relevant implications on the dynamics of Hamiltonian systems. 
Among the most important are: 
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Canonical transformations. A change of variables z = ipiy) is canonical, namely it 



preserve the structure of (1.1), if and only if it is symplectic. Canonical transforma- 



tions were known from Jacobi and used to recast (1.1) in simpler form. 



Volume preservation. The flow 0t of a Hamiltonian system is volume preserving in 
phase space. Recall that if ^ is a (suitable) domain of M^™, we have: 



vol(\/) = / dy, vol(0i(\/)) = / dy 

'V J<f>t{y) 



det^'^'(^) 



dy 



dy. 



However, since ^^^^ = A{t) is a symplectic matrix, from Aiff'" JA{t) = J it follows 
that dei{A{t)Y = 1 for any t and, hence, vol(0t(V)) = vo\{V). 

More in general, volume preservation is a characteristic feature of divergence-free 
vector fields. Recall that the divergence of a vector field / : M" — t- is the trace of 
its Jacobian matrix: 

diyfiy) = ^ + ^ + . . .+ 



dyi dy2 "' dy^' 

The vector field JVH associated with a Hamiltonian system has zero divergence. 
In fact, considering that JVH = [|^, . . . , . . . , -l^l"^ we obtain 

divVH- -0 

dqidpi dqmdpm dpidqi dpmdqm 

since the partial derivatives commute. An important consequence of the previous 
property is Liouville's theorem, which states that the flow (pt associated with a 
divergence-free vector field / : M" — i- M" is volume preserving. 

The above properties and the fact that symplecticity is a characterizing property 
of Hamiltonian systems somehow reinforces the search of symplectic methods for their 
numerical integration. A one-step method 

yi = ^hiyo) 

is per se a transformation of the phase space. Therefore the method is symplectic if is 
a symplectic map. An important consequence of symplecticity in Runge-Kutta methods 
is the conservation of all quadratic first integral of a Hamiltonian system. 



A first integral for system (1.1) is a scalar function I{y) which remains constant if 



evaluated along any solution y{t) of (1.1): I{y{t)) = I{yo) or, equivalently, 

VI{yYJVH{y) = 0, for any y. 

A quadratic first integral takes the form I{y) = y^Cy, with C a symmetric matrix. 

As previously seen, the most noticeable first integral of a Hamiltonian system is the 
Hamiltonian function itself. It is worth noticing that while in the continuous setting en- 
ergy conservation derives from the property of symplecticity of the flow (see, e.g., [55]). 
as sketched above, the same is no longer true in the discrete setting: a symplectic inte- 
grator is not able to yield energy conservation in general. Consequently, devising energy 
conservation methods form an important branch of the geometric integration. 

Symplectic methods can be found in early work of Grobner (see, e.g., [38]). Symplectic 
Runge-Kutta methods have been then studied by Feng Kang [M], Sanz Serna [51], and 
Suris [57]. Such methods are obtained by imposing that the discrete map, associated with 
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Figure 1.1: Level curves for problem (1.7)-(1.9). 



a given numerical method, is symplectic, as is the continuous one. In particular, in 
an easy criterion for simplecticity is provided, for an s-stage Runge-Kutta method with 
tableau given by 



c 


A 







where, as usual, c = (cj 
of the weights, and A = 



e 



I* is the vector of the abscissae, b = G M 
G M*^'* is the corresponding Butcher matrix. 



(1.5) 
is the vector 



Theorem 1.1 (1.5) is symplectic if and only if, by setting B = diag(b), one has: 

BA + A^B = bb^. 



:i.6) 



Since for the continuous map symplecticity implies energy-conservation, though this is no 
more true for the discrete one, then one expects that at least something similar happens 
for the discrete map as well. As a matter of fact, under suitable assumptions, it can be 
proved that, when a symplectic method is used with a constant step-size, the numerical 
solution satisfies a perturbed Hamiltonian problem, thus providing a quasi-conservation 
property over "exponentially long times" Even though this is an interesting feature, 
nonetheless, it constitutes a somewhat weak stability result since, in general, it does not 
extend to infinite intervals. 

Moreover, the perturbed dynamical system could be not "so close" to the original one, 
meaning that, if the stepsize h is not small enough, the perturbed Hamiltonian could not 
correctly approximate the exact one. As an example, consider the problem defined by the 
Hamiltonian [16] 

H{q,p) =p^ + {(3qf + a{q + p)^-. (1.7) 

The corresponding dynamical system has exactly one (marginally stable) equilibrium at 
the origin. Let us select the following parameters 

/3 = 10, a = l, n = 4, (1.8) 



and suppose we are interested in approximating the level curves of the Hamiltonian (shown 
in Figure 1.1) passing from the points 



(9o,Po) 



:i.9) 
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Figure 1.2: 2-stage Gauss method, h = 10 , for approximating problem (1.7)-(1.9) 



This can be done by integrating the trajectories starting at such initial points, for the 
corresponding Hamiltonian system but, if we use the symplectic 2-stage G auss method, 
with stepsize h = 10~^, we obtain the phase portrait depicted in Figure 1.2 which is 
clearly wrongj^ 

A way to get rid of this problem is to directly look for energy- conserving methods, 
able to provide an exact conservation of the Hamiltonian function along the numerical 
trajectory. 

The very first attempts to face this problem were based on projection techniques 
coupled with standard non conservative numerical methods. However, it is well-known 
that this approach suffers from many drawbacks, in that this is usually not enough to 
correctly reproduce the dynamics (see, e.g., |10l p. 111]). 

A completely new approach is represented by discrete gradient methods, which are 
based upon the definition of a discrete counterpart of the gradient operator, so that energy 
conservation of the numerical solution is guaranteed at each step and for any choice of 
the integration step-size [371 EI] ■ 

A different approach is based on the concept of time finite element methods [55] , where 
one finds local Galerkin approximations on each subinterval of a given mesh of size h for 
the equation (1.1). This, in turn, has led to the definition of energy- conserving Runge- 



Kutta methods El ESI I 

A partially related approach is given by discrete line integral methods [IHl SZl HH], 
where the key idea is to exploit the relation between the method itself and the discrete 
line integral, i.e., the discrete counterpart of the line integral in conservative vector fields. 
This tool yields exact conservation for polynomial Hamiltonians of arbitrarily high-degree, 
and results in the class of methods later named Hamiltonian Boundary Value Methods 
(HBVMs), which have been developed in a series of papers [lOlIIIlEllISllIlIISlIIZllIH]. 

Another approach, strictly related to the latter one, is given by the Averaged Vector 
Field method [531 [29] and its generalizations [39], which have been also analysed in the 
framework of B-series [30] (i.e., methods admitting a Taylor expansion with respect to 
the step-size), see e.g., [12]. 



^Additional examples may be found in reference [9]. 
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Further generalizations of HBVMs can be also found in [T^ [201 El El ES]- 



1.2 Discrete line integral methods 



The basic idea which such methods rely on is straightforward. We shall first sketch it in 
the simplest case, as was done in |36], and then the argument will be generalized. Assume 



that, in problem (1.1), the Hamiltonian is a polynomial of degree u. Moreover, starting 
from the initial condition i/q we want to produce a new approximation a.t t = h, say yi, 
such that the Hamiltonian is conserved. By considering the simplest possible path joining 
j/o and yi, i.e., the segment 

a{ch) = cyi + {l-c)yo, ce[0,l], (1.10) 

one obtains: 



H{yr)-H{y,) 



H{a{h))-H{am 
h 

VH{a{t)fa\t)dt 



h [ VH{a{ch)f(x'{ch)dc 
Jo 

h / VH{cyi + (1 - c)yof{yi - yo)dc 
Jo 



-I T 



VH{cyi + (1 - c)?/o)cic 



{yi - yo) = 0, 



provided that 



yi = yo + hJ / VH{cyi + (1 - c)yo)dc. 
Jo 



[1.11) 



In fact, due to the fact that J is skew symmetric, one obtains: 

-1 nT 



h 



-1 



VH{cyi + (1 - c)yo)dc 



uo 



{yi - yo) 



f VH{cy, + {l-c)yo)dc J [ VH{cyi + {l 
Jo J Uo 



c)yo)dc 



0. 



If if G IIj,, then the integrand at the right-hand side in (1.11) has degree v — 1 and, 
therefore, can be exactly computed by using, say, a Newton-Cotes formula based at p 
equidistant abscissae in [0,1]. By setting, hereafter, 

/(■) = JVii(-), (1.12) 

one then obtains 



i=l 



where 



yi = yo + h^ l3if{ciyi + (1 - Ci)yo) = yo + hfiYi) 

i=l 

Yi = (t(q) = c^yl + (1 - Ci)yo, i 



i - 1 



U — 1 

and the are the quadrature weights: 

,1 u 



1,...,!/, 



:i.i3) 
:i.i4) 



n 



t — c 



^dt, 
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Some examples: 

• when u = 2, one obtains the usual trapezoidal rule, 

yi = yo + ^ ifiVo) + fivi)) 

• when u = 3, one obtains the following fomula: 

i/i = yo + ^(/M+4/ + 

• when u = 5, one obtains the formula: 

m = .0+^ (7/fe„) + 32/ (?^) + 12/ (^) + 32/ 



4 



7/(1/1 



The above fromulae were named s-stages trapezoidal rules in [16]. They provide exact 
conservation for polynomial Hamiltonian functions of degree no larger than 2 f or all 



u > 1. Their order of accuracy can be easily determined by recasting (1.13)-(1.14) as a 
z/-stage Runge-Kutta method: 



cb 



T 



with c = (ci, . . . , Cy)"^ and b = {bi, 



(1.15) 



which satisfies some of the usual simplifying assumptions (see, e.g., [HI p. 71]) for an 
s-stage Runge-Kutta method ( see ([1.5[) with coefficients 6j, Cj, Ojj, j — 1, . . . , s: 



q-l ^ 1 



g = l,.. 



g ^ J' 



1, . . . ,S, 



In such a case, in fact, the following result holds true. 

Theorem 1.2 (Butcher, 1964) If a Runge-Kutta method satisfies conditions B{p), C{ri), 
and -D(C)j with 

p < min{?7 + ( + 1, 2(?7 + 1)}, 

then it has order p. 



As a matter of fact, 5(2) and C(l) turn out to be satisfied, for (1.15), thus resulting in 
a second order method. In more details: 

• the quadrature is exact for polynomials of degree 1, so that -8(2) holds true; 

• moreover, by setting e = (1, . . . , 1)^, one has 

cb^e = c ^ C(l). 



Remark 1.1 It is worth mentioning that, even though (1.15) is formally a u-stage im- 



plicit Runge-Kutta method, nevertheless the actual size of the generated discrete problem 
consists of only one nonlinear equation, in the unknown yi, as the above examples clearly 
show. The mono-implicit character of these methods comes from the fact that their coef- 
ficient matrix has rank one. 
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1.3 Generalizing the approach 

The next step is to generalize the above approach, where we have assumed that the path 
a{ch), defined in (1.10), is a hnear function. Now we consider a polynomial path a of 
degree s > 1. Having fixed a suitable basis {Pq, . . . , Pg-i} for Hs^i, one can expand the 
derivative of a as 



s-l 



:i.i6) 



a'{ch) = J2Pjichj^ CG [0,1], 

j=0 

for certain set of coefficients {7^} to be determined. By imposing the initial condition 

0^(0) = yo, 

one then formally obtains 

a{ch)=yo + hy2 i",(x)dx7„ cg[0,1], (1.17) 

with the new approximation given by yi = a{h). Energy conservation may be obtained 
by following a similar computation as before, namely 



H{y^) - H{yo 



H{aih)) - H{a{0)) 

h 

VH{a{t)fa'{t)dt 



h / VH{a{ch)fa\ch)dc 
Jo 

h / Vif((T(c/i))^ VP,(c)7jdc 

-'0 ,_n 



j=0 



j=0 



VH{a{ch))Pj{c)dc 



I3 



J=0, 



:i.i8) 



provided that the unknown coefficients {7^} satisfy 

-fj = r]jJ [ VH{a{ch))Pj{c)dc, 
Jo 

for a suitable set of nonzero scalars 770, ... , r]s-i- The new approximation is then given by 
plugging (1.18) into (1.17): 

s-l „i „i 

y, = aih) = yo + hJ2vj Pjix)dx P,(r)/(a(r/i))dr. (1.19) 

j=0 "'O "'O 



As before, assume H E Hi,. Then, the integrands in (1.18) and (1.19) have at most 
degree [u — l)s + u — 1 = us — 1. Therefore, by fixing a suitable set of k abscissae 
< ci < . . . < Cfc < 1, and corresponding quadrature weights {bi, . . . ,bk}, such that 
the resulting quadrature formula is exact for polynomials of degree us — 1, the integrals 
in (1.18) and (1.19) may be replaced by the corresponding quadrature formulae, which 

k 

7i = Vj hfHc,h))P,{c,), J = 0, . . . , s, (1.20) 



yields 



i=l 
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and 



y, = aih)=yo + hJ2vj / PAx)dxJ2b^PJic^)fiaic,h)), (1.21) 



respectively. By setting, as before, 

Yi = a{cih), i = i,...,k, 



one then obtains: 



Yi = yo + hJ2 



k 



bj'^r]iPi{cj) / Pi{x)dx 



fiYj] 



= yo + h^aijf{Yi), 2 = l,...,/c, 

k 

yi = yo + hY^ 

1=1 



s-i „i 

bi^riiPiici) / P^(x)dx 



:i.22) 



:i.23) 



i=l 



We are then speaking of the following /c-stage Runge-Kutta method: 



A = (aij) e 



sfexfc 



with c = (ci, . . . ,CkY 



h = {h,...,hf, 



with aij, hi defined according to (1.22) and (1.23), respectively. 



:i.24) 



In so doing, energy conservation can always be achieved, provided that the quadrature 
has a suitable high order. For example, we can place the k abscissae {q} at the k 
Gauss-Legendre nodes on [0, 1] thus obtaining maximum order 2k. In such a case, energy 
conservation is guaranteed for polynomial Hamiltonians of degree no larger that 



V < 



2k 



[1.25) 



However, it is quite difficult to discuss the order of accuracy and the properties of 



the fc-stage Runge-Kutta method (3.14), when a generic polynomial basis is considered. 



As matter of fact, different choices of the basis could provide different methods, having 
different orders. As an example, fourth-order energy-conserving Runge-Kutta methods 
were derived in |171 HH], by using the Newton polynomial basis, defined at the abscissae 
{ci}. We shall see that things will greatly simplify, by choosing an orthonormal polynomial 
basis. 

Remark 1.2 It is worth noticing that we can cast in matrix form the Butcher tableau of 



the k-stage Runge-Kutta method (3.14) by introducing the matrices 

Po(ci) ... Ps^lici] 



v.. 



nkxs 



Po(Cfc) ••• Ps-liCk) 
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/£Po{x)dx ... £P,^^{x)dx\ 
\J^'Po{x)dx ... J^' P,_^{x)dx J 




and the row vector 

= { Po{x)dx ... J^' Ps-,{x)dx ) , (1.26) 

as follows: 



c 









which will be further studied later. 



Chapter 2 

Background results 



In this chapter, we state a few prehminary results concerning Legendre polynomials and 
differential equations, for later reference. This chapter is based on references flUl ITi] . 



2.1 Legendre polynomials 

The following polynomials, denoted by Pi, are the Legendre polynomials shifted on the 
interval [0, 1], and scaled in order to be orthonormal: 

degPi = ^, / Pi{x)Pj{x)dx = 6ij, V^,j>0, (2.1) 
Jo 

where 6ij is the Kronecker symbol (its value is 1, when i = j, and 0, otherwise). As any 
family of orthogonal polynomials, they satisfy a 3-terms recurrence, which is given, in this 
specific case, by: 

Po{x) = 1, Pi{x) = V3i2x - 1), 

(2.2) 

P,+i(a;) = {2x - l)^^./55?P,(x) - —J'^^P^_,{x), i > 1. 
^ ^ ' i + 2i + l ^ ' z + 1 V 2z - 1 ^ ^' 

We recall that the roots {ci, . . . ,Cfc} of Pk{x) are all distinct and belong to the interval 
(0, 1). Thus they may be identified via the following conditions: 

-Pfc(ci) = 0, i = l,...,k, with < ci < . . . < Cfc < 1. (2.3) 

It is known that they are symmetrically distributed in the interval [0,1]: 

Ci = l-Ck-i+i, i = l,...,k. (2.4) 

They are referred to as the Gauss-Legendre abscissae on [0, 1], and generate the Gauss- 
Legendre quadrature formula of order 2k, namely an interpolating quadrature formula 
which is exact for polynomials of degree no larger than 2k — 1. In fact, if p{x) G 112^-1, 
then it can be written as 

p{x) = q{x)Pk{x) + r{x), q{x), r{x) G Ilfc.i. 

Consequently, 

.1 .1 

p{x)dx = / [q{x)Pk{x) + r{x)]dx 
Jo 

q{x)Pk{x)dx + / r{x)dx = / r{x)dx. 
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since Pfc(x) is orthogonal to polynomials of degree less than k (see (5.6)). On the other 
hand, for the quadrature formula (q, bi), with the quadrature weights given by 

,1 k 



n 



dx, 



i — 1, . . . , k, 



one obtains: 

k 



i=l 



i=l 



q{ci) Pk{ci) +r{ci) 



r(x)dx, 



due to the fact that any quadrature based at k distinct abscissae is exact for polynomials 
of degree no larger than k — 1. As a matter of fact, for such a quadrature formula, for 
any function / G C^'^([0, 1]), one has 



•^0 i=l 



(2.5) 



for a suitable ( G (0, 1), and with pk independent of /. More in general, if the quadrature 
would have order q < 2k, one would obtain 



1 ^ 
[ /(a;)dx = V6,/(c,) + A,, A, = pfc/(^)(C), 
•^0 i=l 



(2.6) 



with ( and pk defined similarly as above, thus showing that the formula is exact for 
polynomials of degree no larger than q — 1. 

In particular, in the sequel, we shall need to discuss the case where the integrand in 
(2.6) has the following form. 



/(r) = P,(r)/(r/.), rG [0,1], 



(2.7) 



with Pj the jth Legendre polynomial. The following result then holds true. 

Lemma 2.1 Let f G C^'^\ q being the order of the given quadrature formula {ci,bi) over 
the interval [0, 1]. Then 



f P,{T)f{Th)dT - kPj{c,)f{c,h) = 0{h^-^] 

•^0 i=i 



J = 0, . . . , g. 



Proof. The thesis follows from (2.6), by considering that 



d<? 



P,{T)f{Th) ^ [P,{r)f{rh)] 



(9) 



i=0 



E ( ? ) Pf{r)f^-'^{TK)h!'-^ = Oih"-^] 



since Pj^\r) = 0, for i > j.D 



We also need a further result concerning integrals with integrands in the form (2.7), 
which is stated below. 
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Lemma 2.2 Let G : [0, h] — )■ V , with V a suitable vector space, a function which admits 
a Taylor expansion at 0. Then 



[ Pj{T)G{Th)dT = 0{h^), J>0. 

Jo 



Proof. One obtains, by expanding G{Th) at r = 0: 

k>0 ' k>0 

where last but one equality follows from the fact that 



[ P, (r)r'=dr = 0, for k < j.D 
Jo 



2.2 Matrices defined by the Legendre polynomials 

The integrals of the Legendre polynomials are related to the polynomial themselves as 
follows. For all c G [0, 1]: 

^ Po{x)dx = 6^i(c) + ^Po(c), 

/ P,ix)dx = i^+lP^+l{c)-i,P,_^{c), I > I, (2.8) 

Jo 



Remark 2.1 From the orthogonal conditions (5.6), and taking into account that Po{x) 
1, one obtains: 



[ Po{x)dx = l, [ Pj{x)dx = 0, Vj > 1. 
Jo Jo 



(2.9) 



Legendre polynomials possess the following symmetry property: 

P,{c) = {-iyP,il-c), CG [0,1], j>0. (2.10) 
Consequently, their integrals share a similar symmetry: 



[ ' P,{x)dx = i-iy [ ' P,{x)dx, Ci,C2G[0,l], j>0. (2.11) 

J ci J1-C2 

In the sequel, we shall use the following matrices, defined by means of the Legendre 
polynomials evaluated at the k > s abscissae (2.3)|^ 



Po(Cl) ... Ps-M) 

Vs= \ \ ; I G M'^^^ (2.12) 

Po(cfc) ••• Ps-liCk) 



^They have been formally introduced, for a generic polynomial basis, at the end of the previous chapter^ 
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and 



Po(a:)dx ... /;^P.-i(a:)d 



X 



ikxs 



/;'^-Po(x)dx ... /;^p.-i(x)d 



X 



Because of the (2.8), they are related by the following relation: 



Xs — Vs+lXg 



e^-i 



We also set 



tikxk 



... 



(2.13) 



(2.14) 



(2.15) 



the diagonal matrix with the corresponding Gauss-Legendre weights. The following simple 
properties then hold true. 



Theorem 2.1 Matrices (2.12) and (2.13) have full column rank, for all s = l,...,k. 

Moreover, 

VjnVs+i= {Is ). (2.16) 

Proof. By considering any set of s < A; rows of Vs, the resulting sub-matrix is the Gram 
matrix of the s linearly independent polynomials Pq, . . . , Pg-i defined at the corresponding 
s (distinct) abscissae. It is, therefore, nonsingular and, then, Vs has full column rank. 
Moreover, when s = k + 1, one has 



Vk+i = {Vk ) 



(2.17) 



since the last column contains Pk{ci) = 0, i = 1, . . . , k. As a consequence, because of 
(2.14), for matrix one obtains: 

• when s < k, then both Vs+i and Xg have full column rank and so has X^; 



• when s = k, then from (2.17) it follows that 

Xfc = Vk+iXk = VkXk, 

and both Vk and Xk are nonsingular]^ 

Concerning (2.16), one has, by considering that the quadrature formula (cj, bi) is exact for 
polynomials of degree no larger that 2fc — 1 > 2s — 1, and setting ej e and e^- G 
the ith and jth unit vectors: 

k „i 
eJVjnVs+iej = Y^hPi_i{ci)Pj^i{c,) = / P'_i(x)P,_i(x)dx = S^j, 
e=i -^0 

Vz = 1, . . . , s, and j = 1, . . . , s + 1. □ 

From the previous theorem, the following result easily follows. 
Corollary 2.1 When k = s, then V^^ = V'^Vt. 



^Indeed, it can be proved that Xk is nonsingular for all fc > 1. 



2.3. ADDITIONAL PRELIMINARY RESULTS 
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2.3 Additional preliminary results 

In order to deal with the analysis of the methods, we need the following perturbation result 
concerning the solution of the initial value problem for ordinary differential equations: 

y'it) = fiy{t)), t > to, y{to) = yo, (2.18) 

whose solution will be denoted by to, l/o); in order to emphasize the dependence on 
the initial condition, set at (to,l/o)- 

Associated with this problem is the corresponding fundamental matrix, to), satis- 
fying the variational problem (see also ( 1.4[ )) 

$'(t, to) = JMt; to, yo))Ht, to), t > to, $(to, to) = /, 

where the derivative (i.e., ') is with respect to t, and Jf{y) is the Jacobian matrix of f{y). 
The following result then holds true. 



Lemma 2.3 With reference to the solution y{t;tQ,yo) of problem (2.18), one has: 

(i) -^y{t;to,yQ) = ^{t,to); {li) —y{t;to,yo) = -^{t,to)f{yo). 

Proof. Let us consider a perturbation 6yo of the initial condition, and let y{t; to, yo + Syo) 
be the corresponding solution. Consequently, 

y'it;to,yo + 6yo) = f{y{t;to,yo + 6yo)) 

= fijjit; to, yo)) + Jf ivit; to, yo)) [yit; to, yo + 6yo) - y{t; to, yo)] 

^ V ' 

= + O {\\y{t-to,yo + 5yo)-y{t-to,yo)\\') ■ 

Therefore, by setting 

z{t) = yit; to, yo + Syo) - y{t; to, yo), 
one obtains that, at first order (as is the case, when we let 6yo — )■ 0), 

z'it) = Jf{y{t;to,yo))z{t), z{to) = 5yo. 

The solution of this linear problem is easily seen to be 

z{t) = ^t,to)6yo 

and, consequently, 

-7^y{t;to,yo) = -wT^zit) = '^{t,to), 
oyo d{6yo) 

i.e., the part (?) of the thesis. 

Concerning the part (ii), let consider a scalar e ^ and observe that, by setting, 
y{t) = y{t;to,yo), then 

y{t;to + e,yo) = y{t-e). 
Consequently, at first order, the solution of the perturbed problem 

y'it) = f {yit)), t>to + 6, y{to + e)=yo, (2.19) 

coincides with that of the problem 

y'{t) = f{y{t)), t>to, y{to) = yo{e) = yo-ef{yo). (2.20) 
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Letting £ — > 0, one then obtains: 

= -f{yo) 

d d d 

g^y{t;to,yQ) = —y{t;to,yo) g^Voi^) = -*(^, ^o)/(z/o)- 
This concludes the proof. □ 



Chapter 3 

A Framework for HBVMs 



In this chapter, we provide a novel framework for discussing the order and conservation 
properties of HBVMs, based on a local Fourier expansion of the vector field defining 
the dynamical system. In particular, this approach allows us to easily discuss the linear 
stability properties of the methods. The material in this chapter is based on [HI [181 [16]. It 
is worth noticing that an interesting extension of this approach has been recently proposed 
in m. 



3.1 Local Fourier expansion 

Legendre polynomials, previously introduced, constitute an orthonomal basis for the func- 
tions defined on the interval [0, 1]. Therefore, we can formally expand the second member 
of (1.1) over the interval [0, h] as follows (we use the notation ( |1.12[ )): 

/(2/(c/i)) = 5^P,(c)7,(y), CG [0,1], (3.1) 
i>o 

where 

"1 



7,(2/)=/ P,{T)f{y{rh))dT, j > 0. (3.2) 







The expansion ( |3.l[ )-( [3^ ) is known as the Neumann expansion of an analytic function]^ 
and converges uniformly, provided that the function g{c) = f{y{ch)) has continuous 
second derivative for sake of simplicity, hereafter we shall assume g{t) to be analytic. 
In so doing, we are transforming the initial value problem 

y'{t) = f{y{t)), te[0,h], 2/(0) =2/0, (3.3) 

into the equivalent integro- differential problem 



(r)/(2/(r/i))dr, cG [0,1], 2/(0) = 2/0- (3.4) 



In order to obtain a polynomial approximation of degree s to (3.3)-(3.4), we just truncate 
the infinite series to a finite sum. The resulting initial value problem is 

a\ch) = 5^P,(c) / P,{T)f{a{rh))dr ^ P,(c)7,(a), cG [0,1], 
j=0 j-^o 

(t(0) = 2/0, (3.5) 



^E.T. Whittaker, G.N. Watson, A Course in Modern Analysis, Fourth edition, Cambridge University Press, 
1950, page 322. 

^E. Isaacson, H.B. Keller, Analysis of Numerical Methods, Wiley & Sons, 1966, page 206. 
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and its solution evidently defines a polynomial a G 11^. Integrating both sides of (3.5) 
yields the equivalent formulation 



s-l 



a{ch) = yo + h^l Pj{x)dx-fj{a), cG [0,1], 

j=0 



(3.6) 



where the notation (3.2) has been used again. One easily recognizes that (3.5) defines 
the very same expansion (1.16)-(1.18) seen before (with all rij = 1). Consequently, such 
a method is energy-conserving, if we are able to exactly compute the integrals providing 
the coefficients 7j(cr) at the right-hand side in (3.5). From Remark 2.1, one obtains that 

-h 

(3.7) 



^(h) =yo+ /(a(r/i))dr. 

^0 

Let now discuss the order of the approximation cr[h) ^ y{h). 



Lemma 3.1 Let 7j(cr) be defined according to (3.2). Then 7^(0") = 0{h^). 



Proof. The proof follows immediately from (3.2), by vrtue of Lemma 2.2[ □ 

We are now able to prove the following result. 
Theorem 3.1 a{h) - y{h) = 0{h'^'+^). 



Proof. Denoting by y{t;tQ,yo) the solution of problem (2.18) and considering that o"(0) 
yo, by virtue of Lemmas 2^ and 3J^ one has: 

a{h)-y{h) = y{h; h, a{h)) - y{h;0,yo) 

= y{h;h,aih))-yih;0,aiO)) = 



dt 



y{h-t,a{t))dt 



^ y{h-e,a{t)) + -^y{h;t,oj) ,a'it)]dt 



89' 



dco' 



U]=a(t) 



[-$(/i,t)/(a(t)) + $(t,to)^'(t)] dt 
^{h,t)[-f{a{t)) + a\t)]dt 

h [ <!>{h,Th)[-f{(r{Th)) + a\Th)]dT 
Jo 



-h / (!>{h,rh) 



s-l 



5^P,(r)7,(a)-5^P,(r)7,(a) 



.i>o 



j=0 



dr 



-h [\ih,Th)Y,PjirhA^)dT 



j>s 

= G{Th) 



^{h,Th) Pj(r)dr 



= o(y) 



= 0(ho) 



j>s 



3.1. LOCAL FOURIER EXPANSION 
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We observe that, unless we can compute exactly the integrals defining the {7^(0")} in (3.5) 
(which is the case, for example, when / is a polynomial or in very special situations), 
(3.5) is not yet an operative method, but rather a formula. In order to obtain a numerical 
approximation procedure, we need to approximate those integrals by means of a suitable 
quadrature formula, which we define at A; > s Gauss abscissae in [0, 1] defined in (2.3). As 
was seen in Section 2.1, the corresponding quadrature formula (cj,6j) has order q = 2k, 



that is it is exact for polynomials of degree no larger that 2k — 1. By recalling Lemma 2.1 



let then approximate the integrals in (3.5) by means of a quadrature (cj, 6,) over k distinct 
abscissae. Consequently, in place of a defined by (3.5) or (3.6), we shall compute the new 
polynomial m G such that: 



s-l 



u\ch) = J2P,{c)Y,biPj{ci)f{u{ceh)), cG[0,l], 
j=0 i=\ 
m(0) = yo, 



that is, 



u{ch) =yQ + hJ2 / Pj{x)dx J2 biPj{ci)f{u{ceh)), c G [0, 1], 
with the new approximation given by: 

k 

= u{h) = yo + h^bif{u{cih)). 



(3.8) 



(3.9) 



(3.10) 



4 = 1 



If the quadrature formula {ci,bi) has order g, then, by virtue of Lemma 2.1 and taking 
into account (3.2), one obtains 

1 ^ 
7,(n) = [ P,{T)f\u{Th))dr = J2biPAce)fiuiceh)) - A,(/i), (3.11) 

A,(/i) = Oih^'^), J = 0,...,q. 



Consequently, we can rewrite the first equation in (3.8) in the following equivalent form: 

s-l 



U 



'{ch) = Y,PA^)[iA^)-^m^ cG[o,i]. 



(3.12) 



This allows us to derive the following result, which we state for a generic quadrature of 
order q. 

Theorem 3.2 yi — y{h) = 0{W^^), where p = min{g, 2s}. 



Proof. The proof proceeds on the same line as that of Theorem 3.1 



yi-y{h) = u{h)-y{h) = y{h;h,u{h)) - y{h;0,u{0)) 



dt 



y{h-t,u{t))dt 



' ^ y{h-9,u{t)) +^y{h-t,u) u'{t)]dt 

e=t ou ui=u{t) 

1 



39' 



<!>{h,t)[-f{u{t))+u'{t)]dt = h / <!?{h,Th)[-f{u{Th))+u'{Th)]d 



IT 
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-h [ (!>{h,Th) 
Jo 

h [\ih,Th)Y,P,ir)A,iu)dT-h l\{Krh)Y,PAr)lM)^r 

Jo .-n Jo 



s-1 
j=0 



j=0 

= G{Th) 



^{h,Th) Pj{T)dT 



= 0{hi-3) 

j>s 



= G{Th) 



= 0{h3) 

7j(M) 



— V 

--O(hi) 



--O{ho) 



0{h''+^) + h^Oih^^) = 0(/iP+^), p = mm{q,2s}.D 

j>s 



Definition 3.1 The method ( 3. 12}-(3.10} defines a Hamiltonian Boundary Value Method 
(HBVM) with k stages and degree s, in short HBVM{k, s). 



As a consequence, by setting the abscissae at the k Gauss points (2.3), the following result 
holds true. 



Corollary 3.1 By choosing the k abscissae {cj} as in (2.3), a HBVM{k, s) method has 
order 2s, for all k > s. 



3.2 Runge-Kutta form of HBVM(A;, s) 

Before studying the conservation properties of the methods, let us derive the Runge-Kutta 
formulation of HBVM(fc, s). The basic fact is that, at the right-hand sides of equations 
(3.9)-(3.10), one only requires to know the value of the polynomial u at the abscissae 

i = 1, 



{cih}. Consequently, by setting 

Yi = u{cih), 
one obtains: 



Y 



yi 



yo 



yo 



i=i 

k 



s-1 



1=0 



+ hY,a.,f{Y,), 



1, . . . , /c. 



k 



yo 



hY^hf{Y,). 



1=1 



In other words, we have defined the following fc-stage Runge-Kutta method: 



with (see (3.13)), 

c = (ci, . . . ,Cfc)^, b = (bi, 



and A = (aij) G 



nkxk 



(3.13) 



(3.14) 



The Butcher tableau (3.14) defines the Runge-Kutta shape of a HBVM(k, s) method. We 
can easily derive a more compact form for the Butcher array A in (3.14). 



3.2. RUNGE-KUTTA FORM OF HBVM{K, S) 



21 



Theorem 3.3 A = XsVjVt, with the matrices Xg^Vs,^ defined according to (2.12)-(2.15) 

the z-th and j-th unit vectors, one obtains: 



Proof. By setting ej, e, G 



( /;^Po(a:)dx 

s-l 



£=0 



6, 



according to (3.13). □ 



Consequently, the Buther tableau (3.14) can be casted as: 



c 









or, equivalently, by taking into account (|2.14|), 

c 



(3.15) 



(3.16) 



Remark 3.1 We observe that the Rung e-Kutta form (3.15) of a HBVM{k, s) method is 



simplified, with respect to that sketched in Remark L2 for a discrete line-integral methods 
defined by using a general polynomial basis. In particular, the diagonal matrix Ag is now 
authomatically fixed, in order to maximize the accuracy of the method, and the vector 
of the quadrature coincide with that used for approximating the integrals involved in the 
coefficients of the polynomial u. 



3.2.1 HBVM(s,s) 

In the case k = s, the matrices X,,Po,n e W^^. 



Moreover, the following results follows 



immediately from Theorem 2.1 and Corollary 2.1 



Consequently, in such a case, we can write the Butcher tableau (3.16) as that of the 
following s-stage method, 

c 



(3.17) 

which is the W -transform defining the s-stage Gauss-Legendre Runge-Kutta collocation 
method [HI p. 79], which has order 2s. In this sense, in the case k > s, HBVM(A;,s) 
can be ragarded as low-rank generalizations of the s-stage Gauss method. Indeed, the 
following result holds true. 

Theorem 3.4 For all k > s the rank of the matrix A = Vs+iXgVjfl is s. Moreover, the 
nonzero eigenvalues coincides with those of the basic s-stage Gauss method. 

Proof. The rank of the matrix Vg+i is s or s + 1 (when k > s), whereas that of matrices 
Xg^Vs is s, and Vt is nonsingular. Therefore, the rank of A cannot exceed s. Moreover, 
from Theorem 2J^, one has 

v'^^AVs = rJnVg+iXgVjnVg = (/, o)x,j, = g w% 
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which is known to be nonsingular. Consequently, rank(y4) = s. Moreover, 

vJnA = vJnVs+iXsVjn = (/, o)XsVjn = XsVjn. 

This means that the columns of QVs span an s-dimensional left invariant subspace of 
A. Therefore, the eigenvalues of Xg will coincide with the nonzero eigenvalues of A. On 



the other hand, from (3.17) one obtains immediately that the eigenvalues of Xs are the 



eigenvalues of the Butcher matrix of the s-stage Gauss method. □ 

This property has been named isospectrality of HBVMs, in [17j. It will be used for the 
efficient implementation of HBVM(/c, s) methods. 



3.3 Energy conservation 



We now consider the issue of energy conservation for HBVM(/c,s) methods. From (3.8) 
(3.10) with / = JVH, we obtain: 



H{yi)-H{yo) = H{uih)) ^ H{u{0)) = / VH{u{t)fu'{t)dt 

Jo 

= h [ VH{u{Th)fu'{Th)dT 

Jo 

„1 s— 1 k 

= h S7H{u{Th)f^Pj{T)^hiPj{ci)JWH{c,h)dT 

.i=0 i=l 



s-1 
j=0 



Pj{T)JVH{u{Th))dT 



— V 

--0{hi) 



J 



Y,biPj{ci)JVH{cih) 



i=l 



= Eh 

Now, two possibilities may occur: 



/ Pj{T)JVH{u{Th))dT = \biPj{ci)JVH{cih): in such case. Eh = 0, so that 

i=l 

energy is exactly conserved. This is the case of a polynomial Hamiltonian of degree 
u no larger that 2k /s; 

nl k 

/ P,j{T)JV H{u{Th))dT = 'S^hiPj{ci)JV H{cih) — Aj{h) : in such a case, by taking 

-'0 i=l 



into account (3.11 ), one obtains that Eh = 0{h'^^'^^), provided that the Hamiltonian 



is suitably regular, as we have assumed. 
We have then proved the following result. 

Theorem 3.5 HBVM{k, s) is energy-conserving for all polynomial Hamiltonian of degree 

2k 

z/<— . (3.18) 



In any other case, H{yi) — H{yo) = Oiji^^^^), even though the method has order s. 



3.4. SYMMETRY 
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Remark 3.2 We observe that: 

• for polynomial Hamiltonians, energy conservation can he always obtained, by choos- 



ing k large enough, by virtue of (3.18); 



even in the case of non polynomial Hamiltonians, energy conservation can be prac- 
tically gained by choosing k large enough, provided that \Eh\, which is 0{h^''^^), is 
within roundoff errors. 



As an example, in Figure we plotted the level curves passing at 



(?o,Po) 



U, -I) 



for the Hamiltonian problem with Hamiltonian 

H{q,p)=p'' + {Pqf + a{q + pf\ 



with parameters: 



/3 = 10, 



a 



n 



4. 



(3.19) 



(3.20) 



(3.21) 



By using the 2-stage Gauss method (fourth-order), with step size h = ^, the obtained 
phase portrait is wrong, as is shown in Figure |3.1[ due to the error in the numerical Hamil- 



tonian, which is shown in Figure |3.2[ Indeed, even though no drift in the Hamiltonian 
occurs, nevertheless it is not negligible, for the problem at hand. 

However, if we use HBVM(3, 2) with the same step-size, the error in the Hamiltonian 
is of sixth-order: this is enough to have a smaller error in the numerical Hamiltonian, as 
is shown in Figure |3^ resulting in a correct phase portrait, as is shown in Figure 3.3 



At last, by using HBVM(8,2) with the same step size, the Hamiltonian error is of the 



order of roundoff errors, as is shown in Figure 3.6, thus allowing a perfect reconstruction 



of the phase portrait, depicted in Figure 3.5 Indeed, since the Hamiltonian (3.20) has 



degree eight, the quadrature is exact, in this case, according to (3.18). 



For sake of completeness, in Fugure 3/7 we also plot the mean error in the numerical 
Hamiltonian, for HBVM(A;, 2) methods, used with the stepsize h = 10~^, for = 2, . . . , 8. 
As one can see, for the largest values of k the error is essentially due to roundoff. 



3.4 Symmetry 



We here prove that, provided that the abscissae {cj} are symmetrically distributed in 
the interval [0,1], as is the case of the Gauss-Legendre nodes (see (2.4)), a HBVM(fc, s) 
method is symmetric. In more detail, if applied to the initial value problem 

y' = f{y), y{o)=yo, 

it provides the approximation yi ^ y{h), then it will provide the same discrete solution, 
as well as the same internal stages, though in reversed order, if applied to the initial value 
problem 

z' = -f{z), z{0) = y,. (3.22) 
For proving this, let us define the following matrices: 



1 \ 



k,k + l,k + 2, 
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Figure 3.4: Hamiltonian error for problem (3.19)-(3.21 ), HBVM(3,2) method, h = 10 ^. 
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Figure 3.6: Hamiltonian error for problem (3.19)-(3.21 ), HBVM(8,2) method, h = W ^. 
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10 
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10 



< 10 

o 



10 



10 



10 



7 



Figure 3.7: Mean Hamiltonian error for problem (3.19)-(3.21 ), HBVM(A;,2) method, k 
2, . . . , 8, by using a stepsize h = 10~^. 



/ 1 



L 



-1 1 



v 



-1 1/ 



/ 1 



nk+lxk+l 



D 



\ 



and, by recalling the vector 1} defined at (1.26), 



Is 

X] 



X, = 1 -1 I G 



)fc+lXf 



Moreover, by setting 



= Co < Ci < ■ ■ ■ < Cfe < Cfe+l = 1, 



(3.23) 



we need to define matrix 

LX. = AX, 



Pj_i(x)dx 



Ci-l 



i = 1, . . . , A: + 1 
j = 1, . . . ,s 



The following properties then hold true, provided that the abscissae are symmetrically 



distributed in the interval [0,1], i.e., by taking into account (3.23), q = 1 — Ck-i, i 
0,...,fc + l: 



JkVs = VsD; 



where the last two properties follow from (2.11) and (2.10), respectively. The discrete 



solution generated by a HBVM(/c, s) method can then be cast in vector form as 

( -e 4+1 ) ® / r = hisV^n ® / /(r). 
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where e G M'^^^ is the unit vector, and 

Left multiphcation by L ® / then gives 

A^IY = hB®I f{Y), 

with 



(3.24) 



-1 1 



A 



B= {0 AlsPjn ) G 



Dfc+lxfe+2 



-1 1 



Since one easily reahzes that 



Jk+lAJk+2 — —A, 

the method would be symmetric provided that 

Jk+lBJk+2 = B. 

In fact, by observing that 

yi 

Z = Jk+2 ®/r=| Jk(^IY 

yo 



is the reversed-time discrete solution, left multiplication of (3.24) by Jk+i ® I then gives: 



= Jk+iA®IY -hJk+iB®I f{Y) 

= Ju+iAJl^^®IY -h.h+iBJl^^®IfiY) 
= -A®IZ -hB®I f{Z). 

That is, the reversed-time vector satisfies the equation 

A®IZ = -hB^I f{Z), 



which consists in applying the HBVM(A;, s) method to problem (3.22), thus providing the 
approximation zi = yo, and using stages Z = Jk ® lY . As matter of fact, one has: 

Jk+iBJk+2 = ( Jk+iAIgV'^QJk ) 

= ( AIsDVjJk^ ) 

= ( AIsD{JkVsfn ) 

= ( AI.D'^VjQ ) = 5, 

and the symmetry of the method follows. 
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3.5 Linear stability analysis 



We now consider the linear stability analysis of HBVM(/c, s): indeed, such methods can be 
defined independently from the problem of energy conservation, by considering a general 
function / in (3.8). As matter of fact, we have seen, in Section 3.2.1, that HBVM(fc, s) 
methods, with k > s can be regarded as a low-rank generalization of the basic s-stage 
Gauss-Legendre method. 

Then, let us apply one such a method to the celebrated test equation 



Setting 



y' = Xy, yiO) = yo 7^ 0, 3ft(A) < 0. 
X = a + y = Xi + ix2, 



the test equation becomes: 



Xi 

X2 

Defining the scalar function 



a —f3 

(3 a 



y(x) 



Xi 
X2 



- T 
-X X 



Ax, x(0) = xo ^ 0. 
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(3.25) 



(3.26) 



the application of a HBVM(A;, s) method for solving (3.25) defines the polynomial a such 
that cr(0) = xq and, moreover 



a'{ch) 



s—l k s—l „\ 

j=0 i=l j=0 

^E^.W f PAr)^V{a{Th))dT. 
Jo 



'T)Aa{rh)dT 



(3.27) 



By considering that cr(0) = xq, and the new approximation is defined by 

xi = a{h), 

one obtains: 

AV(xo) = r(xi)-r(xo) = V{a{h))-V{cr{0)) 



dt 

1 



V{a{t))dt 



VV{a{t)ya'{t)dt 



h fvV{a{rh)rAY,PA^)\fpAoWV{ 
Jo L^o 



a{ch))dc 



dr 



s-l 

ah 

3=0 

s-l 

ah 

j=0 



Pjir)W{a{Th))dT 



Pj{T)a{Th)dT 
Moreover, the following result holds true. 



ahT^. 



(3.28) 
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Lemma 3.2 = 

Proof. Indeed, one has: 

= 



xo = 0. 



a'{ch) = and 



Po{ch) cr(c/i)dc 



[ cr{ch)dc = 0. 
Jo 



From the first equahty one obtains a{ch) = xq and, therefore, from the second equality 
one derives xq = 0. □ 



From (3.28) and Lemma 3.2, the following result then easily follows. 



Theorem 3.6 For all k > s, and for any choice of the nodes, HBVM{k, s) is perfectly 
A-stable, i.e., its stability region coincides with the negative-real complex plane, C~ 



Proof. From (3.28) and Lemma 3.2, one has, by considerdering (lyapV) and that a = 3f?(A): 



||xi||2 = llxolla + a/iP^ < IIX0II2 ^ 3?(A)<0. 

Consequently, a HBVM(fc, s) method turns out to be perfectly A-stable, since its absolute 
stability region coincides with C~, for all A; > s > 1. □ 



Chapter 4 

Implementation of the methods 



In this chapter, we discuss the efficient implementation of HBVM(/c, s) methods. In 
particular, it is clearly shown that their computational cost depends essentially on s, in 
the sense that, for all k > s, the discrete problem turns out to have always block-dimension 
s. A nonlinear iteration procedure, based on the blended implementation of the methods 
is also sketched. The material in this chapter is based on |i8| flOt IT2t ri6[ lll2T| [22| l2i[ 123] . 

4.1 Fundamental and silent stages 



From (3.15)-(3.16), we see that a HBVM(/c, s) method, with A; > s, is defined by a Butcher 
matrix of rank s. Consequently, — s of the stages of the method can be expressed as 
a linear combination of the remaining s stages: we shall, therefore, name fundamental 
stages the latter ones, and silent stages the former ones. For this purpose, let us partition 
the stage vector Y as 

y(2) 



Y 




where, by supposing for sake of brevity that the fundamental satges are the first s-onesQ 

/ n+i 
= : 

V n 

Similarly, we partition matrices Is and Vs, respectively, as 

containing the corresponding rows as those of Y^^^ and Y^'^\ respectively. Moreover, we 
also consider the partition 



^2 



Consequently, by setting e^^^ end e*^^) the unit vectors of length s and k — s, respectively, 
one obtains: 

F« = eW ^yo + hX^^VjU ® / ( ^[^Jj;] ) , (4.1) 



^Indeed, this can be always achieved, by using a suitable premutation of the abscissae. 
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(4.2) 



From (4.1), one then obtains that 

Pin 



= (/.xi^))-^[rW-e(^)®yo], 

which substituted into ( |4.2[ ) gives: 



Consequently, we can rewrite (4.1)-(4.2) as: 



(4.3) 



involving only the fundamental stages, thus confirming that the actual discrete problem, 
to be solved at each time step, amounts to a set of s (generally) nonlinear equations, each 
having the same size as that of the continuous problem. For solving such a problem, one 
could use, e.g., a fixed-point iteration, 



or, if the simplified- Newton iteration. In more details, setting 

(Pp))^fi, ® / / (a® yo + Xf (X«)-V^))] , 



(4.4) 



one then solves, 

[I-hC^Jo] A, = -F(r/')), 

where Jo = Jfivo) and matrix C is defined as follows: 

c = ja) [(pa))^fi, + (Pp))^fiixf )(xi^))-^] 

The following result holds true. 



0,1,..., 



(4.5) 



(4.6) 



Theorem 4.1 The eigenvalues of matrix C, as defined in {4-6), coincide with those of 



matrix Xs defined in (2.1^), that is the eigenvalues of the Butcher matrix of the s-stage 
Gauss method. 



4.2. ALTERNATIVE FORMULATION OF THE DISCRETE PROBLEM 
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Proof. One has: 

= [Is 0]Xs = Xs.D 

Consequently, matrix C has always the same spectrum, independently of the choice of the 
fundamental and silent abscissae^ This, in turn, coincides with the nonzero eigenvalues 
of the corresponding Butcher array (see Theorem 3.4). Nevertheless, its condition number 



is greatly affected from this choice. Clearly, a badly conditioned matrix C would affect 



the convergence of both the iterations (4.4) and (4.5). As an example, in Figures 4.1 and 



4.2 we plot the condition number of matrix C corresponding to the following choices of 
the fundamental abscissae, in the case k > s = 3: 



the first s abscissae of the k ones (Figure 4.1); 



s almost evenly spaced abscissae among the k ones (Figure 4.2 ) 



As one may see, in the first case n{C) grows exponentially with k, whereas it is uniformly 
bounded in the second case. Because of this reason, we shall consider a more favorable 
formulation of the discrete problem itself, which will be independent of the choice of the 
fundamental abscissae. 

4.2 Alternative formulation of the discrete problem 

In order to overcome the previous drawback, the basic idea is to reformulate the discrete 
problem by considering as unknowns the coefficients, say 



7j = ^bePj{ce)f{u{ceh)), 



j = 0, . . . , s - 1, 



of the polynomial approximation defining a HBVM(fc,s) method (see (3.9)). In more 
details, recalling that 



s-l 



Yi = u{cih) = yQ + h'^^j / Pj{x)dx, 



j=Q 



i — 1, . . . , k, 



one may cast the disctere problem as follows: 



7 



7o 



7.-1 



Vjn ®If{e®ya + hXs® Ij) 



with the new approximation given by 



yi = yo + h%. 



(4.7) 



^I.e., the abscissae corresponding to the fundamental and silent stages, respectively. 
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, X 10 




Figure 4.1: Condition number of matrix (4.6), fundamental abscissae fixed at the first s ones. 




Figure 4.2: Condition number of matrix (4.6), fundamental abscissae almost evenly spaced. 
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We observe that (4.7) has always (block) dimension s, whatever is the value of k consid- 
ered. For solving such a problem, one can still use a fixed-point iteration, 

7^+1 = Vjn f {e^yo + hls® I Y) , £ = 0,1,..., 

whose implementation is straightforward. One can also consider a simplified-Newton 
iteration. Setting 

F(7) = 7 - VjQ ^If{e^yo + hIs® I^) , (4.8) 
and, as before, Jq = Jf{yo), it takes the form 

[I-hC®Jo]A' = -F{Y), 7'+^ = y + A^ £ = 0,1,..., (4.9) 

where matrix C is now defined as follows: 

c = vjni, = vJnVs+iXs = (/. o)x, = x,. (4.io) 

Consequently, the iteration ( |4.9[ ) becomes: 

[/-/iX,® Jo]A^ = -F(y), y+i = y + A^ £ = 0,1,.... (4.11) 



Remark 4.1 It is worth noticing that (4-11) holds independently of the choice of the 
k abscissae {q}, the only requirement being the order 2s of the quadrature, so that the 
property Vj^QVs+iXs = (/« 0) holds true. 



Remark 4.2 We observe that both matrices (4-6) CLnd (4-10) share the same eigenvalues 
which, in turn, are the nonzero eigenvalues of the Butcher array of the given HBVM{k, s) 
method (see Theorem 3.4)- 

We observe that, remarkably enough, at each step of the simplified-Newton iteration we 
have to solve a linear system of dimension sm x sm of the form 



[/ - hXs ® Jo] X = r]. 



(4.12) 



whose coefficient matrix is thus independent of k and of the choice of the abscissae. Its 
cost is then approximately given by 

2 

-(sm) flops, 

due to the cost of the corresponding LU factorization. In the next section, we shall 
consider an alternative, iterative, procedure for solving (4.12), able to reduce the cost for 
the factrization to 



-m 



flops. 



4.3 Blended HBVMs 



We now introduce an iterative procedure for solving ( |4.12[ ), which has been already suc- 
cesfully implemented in the computational codes BiM [22j and BiMD [24j| for the numerical 
solution of stiff ODE-IVPs and linearly implicit DAEs up to order 3. 

For this iterative procedure a linear analysis of convergence is provided. To this 
purpose, let us consider the "usual" test equation. 



y' = Xy, 3fJ(A) < 0. 



(4.13) 
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In such a case, by setting as usual q = hX, problem (4.12) becomes the linear system, of 
dimension s, 

(/-gX,)x = r7. (4.14) 

The solution of this linear system is not affected by left-multiplication by CX~^, where 
C > is a free parameter to be chosen later. Thus, we obtain the following equivalent 
formulation of (4.14): 

aX:'-qI)^ = (X-'ri^ri^. (4.15) 
Let us define the weighting function 

9iq)=I{l-Cq)-'. (4.16) 
It satisfies the following properties: 

• 9{q) is well defined for all q G C~, since C > 0; 

• ^(0) = /; 

• ^(9) — O, as g — 7- 00. 



We can derive a further equivalent formulation of problem (4.14), as the blending, with 
weights 6{q) and / — 6{q) of the two equivalent formulations (4.14) and (4.15), thus 
obtaining 

M(g)x = 77(g), (4.17) 



with: 



M{q) 



e{q){I - qX,) + C(/ - e{q)){X:' - ql), 
9iq)v + ai-0iq))X;% 



(4.18) 



Equations (4. 17)- (4. 18) define the blended formulation of the original problem (4.14). The 
next step is now to devise an iterative procedure, defined by a suitable splitting, for solving 
(4.17). To this end we observe that, due to the properties of the weighting function 9{q) 
defined in (4.16), one has: 



M{q) 
M{q) 



q ^ 0, 



\q\ > 1. 



Consequently, N{q) = /(I — (q) ^ M{q), both for q ^ 0, and |g| ^ 1. It is then natural 

ve procedure, for solving 

(iV(g)-M(g))x, + 77(g), 



to define the following iterative procedure, for solving (4.17) 

iV(g)x.+i = 

That is, observing that A^(g)^"'^ 

x,+i = (/ - ^^(g)M(g))x, + ^(g)r7(g), r = 0, 1 



0,1,.... 



(4.19) 



Equation (4.19) defines the blended iteration associated with the blended formulation 
(4.17) of the problem. By considering that the solution, say x*, of (4.17) satisfies also 
( 4.19 ), by setting 

the error at the rth iteration, one then obtains the error equation 

e,+i = (/ - ^(g)M(g))e„ r = 0, 1, . . . . 
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Consequently, the iteration (4.19) will converge to the solution x* of the problem iff the 

(4.20) 



spectral radius of the iteration matrix, 

p(g) = max 1^1, with 

Cea(Z(g)) 



Z{q) = I-e{q)M{q), 



is less than 1, where a{-) denotes the spectrum of the matrix in argument. The set 

B = {qeC : p(g) < 1} 



is the region of convergence of the iteration (4.19). The iteration will be said to be: 

• A-convergent if C D; 

• L-convergent if, in addition, p(g) — )• 0, as g — i- cxd. 

Remark 4.3 A- convergent iterations are then appropriate when the underlying method 
is A-stable. Similarly, L-convergent iterations are appropriate in the case of L-stable 
methods. 

We observe that: 

• Z{0) =0 p(0) = 0; 

• Z{q) — )• O =^ p{q) — )■ 0, as g — oo; 

• Z{q) is well-defined for all q G C~, since ( > 0. 



Consequently, for the blended iteration (4.19) A- convergence and L-convergence are equiv- 
alent to each other. From the maximum modulus theorem, in turn, it follows that this is 
equivalent to requiring that the maximum amplification factor, 

p* = sup p(g) = supp(ix), 

5R{g)=0 xm 

satisfies 

P* < 1- 

For the blended iteration, due to the fact that p(g) — )■ 0, as g - 
Xs is real, so that p(g) = p{q), one has actually to prove that 



p* = maxp(«x) < 1. 

a:>0 



oo, and since the matrix 
(4.21) 



We shall choose the free positive parameter (, in order to minimize p*, so that (4.21) 
turns out to be fulfilled for all s > 1 (see (4.14) and (2.14)). The following result holds 
true. 

Theorem 4.2 p e cr(X,) ^ ^[{Tg^^)^^ ^ ^(^(?))- 



Proof. From ( |4.20[ ), ( |4.18| ), and ( |4.16[ ), one obtains: 

Z{q) = I-e{q)M{q) 

= I - e{q) [e{q){I - gX,) + C(/ - e{q)){X;' - ql)] 

= I - e{qf [(/ - gX,) - eq{X:' " ql)] 

= 0{qf [(1 + Cq^ - Kq)I -I + qXs + CqX;' - (Vl] 

= qe{qfX:' [Xl - 2CX, + C'l] 

= q9{qfX;\X,-CI? 

^ q{X,-Clf[X,{l-Cq)Hy\ 
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from which the thesis easily follows. □ 

As a consequence, one obtains the following result. 

Corollary 4.1 The maximum amplification factor |^ . ) of the blended iteration (4-19) 
is given by: 

p = max — - — — . 

Proof. One has: 

x\i2-C\'^ X |p-CP 
p = max max - — — -- = max max 



The thesis then follows immediately, by considering that 

X 1 



max 



x>o 1 + C2a;2 2C' 
which is obtained at x = C,~^. □ 

We are now in the position to choose the positive parameter C, in order for p* to be 
minimized. This clearly will depend on the eigenvalues of matrix Xs- Since this matrix 
is real, the complex ones occur as complex-conjugate pairs. Consequently, if we set 

Pj ~ Ipil^ "^^5 j ~ ^1 ■ ■ ■ 1 ^1 
we can sort them by decreasing arguments: 

- > 01 > 02 > ... > 0« > - 

due to the fact that 

3f?(/ij)>0, j = l,...,s. 
Moreover, we can neglect the complex conjugate ones, thus obtaining: 

|>0i>...0.>o, ^= 

In addition to this, it turns out that the eigenvalues of matrix Xg also satisfy: 

< \pi\ < ... < \pe\, 



as is shown in Figures |4.3| and |4.4[ in the cases s = 6 and s = 7, respectively. In such a 
case, the following result holds true. 

Theorem 4.3 p* is minimized by choosing 



C = |/ii| = min \p\, (4.22) 



resulting in 



In such a case, one obtains: 



p. 1 Ipi - CI 



(4.23) 

C=lml 



2C l/iil 

= l-cos0i<l. (4.24) 



4.4. ACTUAL BLENDED IMPLEMENTATION 

Table 4.1: Bended iteration of HBVM(A;, s) methods. 
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s 


C 


P* 


P 


2 


0.2887 


0.1340 


0.0774 


3 


0.1967 


0.2765 


0.1088 


4 


0.1475 


0.3793 


0.1119 


5 


0.1173 


0.4544 


0.1066 


6 


0.0971 


0.5114 


0.0993 


7 


0.0827 


0.5561 


0.0919 



Proof. For ( |4.22| )-( [423| ), see pi]. Concerning ( |4.24[ ), one has: 

P* - 



1 - lyUllP _ 


1/^1 


2[(l-cos0i)2 + (sin0i)2] 


2|/^i| \Pl\ 






1 + (cos0i)^ + (sint; 




- 2 cos 01 2 — 2 cos 01 


2 




2 



= 1 — COS 01 . □ 

Consequently, the blended implementation of HBVM(A;, s) methods is always A-convergent 
and, therefore, L-convergent. We can also characterize the speed of convergence when 
g ~ 0, by considering that, from Theorem 4^, Corollary 4A, and Theorem 4^, it follows 
that 



Piq)- 

where the parameter 





Pi - 


Pi\ 


2 


Pi - 




Pl\ 


|2 


\pl\ 


|1 - q\fii 


2 


Pi 





\q\ + 0(|g|^) ^ p\q\ 



P 



\pl - 




Pl\ 


|2 




Pi 





is called the non-stiff amplification factor. In Table |4.1| we list the relevant information 
for the iteration of HBVM(/c, s) methods. 



4.4 Actual blended implementation 

Let us now sketch the blended implementation of HBVMs, when applied to a general, 
nonlinear system, also analyzing its complexity. In the case of the initial value problem 

y=f{y). 1/(0) = G M'", 

the previous aguments can be generalized in a straightworfard way, by considering that 
now the weighting function becomes 

9 = 1,® VL = I^- hCJo, 

where h is the stepsize, ( is the optimal parameter specified in the second column in 



Table 4.1, and Jq is the Jacobian of / evaluated at yo (clearly, we are speaking about the 



first step in the numerical integration). 



From (4.8) and (4.11), we have to solve the onier-mner iteration described in Table 4.2 
Let us analyze its computational complexity (let m denote the dimension of the continuous 
problem and e G be the unit vector), by considering, for each item, only the leading 
term in the complexity and denoting, as 1 flop, an elementary (binary) algebraic /boating- 
point operation. One obtains: 
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0.1 p 
0.08 - 
0.06 - 
0.04 - 
0.02 

' 
-0.02 - 
-0.04 - 
-0.06 - 
-0.08 - 
-0.1 







0.02 



0.04 



0.06 
5i(M,) 



0.08 



0.1 



0.12 



Figure 4.3: Eigenvalues of matrix Xq. 
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0.02 
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o 


o 



0.02 



0.04 
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5i(M,) 
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0.1 
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Figure 4.4: Eigenvalues of matrix X7. 
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• fi: 1 Jacobian evaluation; 

• 9: |m'^ flops for the LU factorization of fl; 

• y^: 2ksm flops; 

• k function evaluations; 

• 1]^: 2ksm flops; 

• z^'^: 2sm? flops; 

• u^''': 2s^m flops; 

• w^'"^: 2s^m flops; 

• A^'^'+i; 4sm2 flops; 

• 7^"^-'^: sm flops. 

Consequenly, this algorithm has a fixed computational cost of 1 Jacobian evaluation and 
|m'^ flops, plus, assuming that v inner iterations are performed, a cost of k function 
evaluations and Aksm + i/(6sm^ + 4s^m) + sm flops per outer iteration. 

A simplified (and sometimes more efficient) procedure is that of performing a nonlinear 
iteration, obtained by performing exactly 1 inner iteration (i.e., that with r = 0) in the 
above procedure, thus obtaining the algorithm depicted in Table 4.3 In such a case, the 
resulting computational cost is obtained as follows: 

• Q: 1 Jacobian evaluation; 

• 6: flops for the LU factorization of Q; 

• y^: 2ksm flops; 

• k function evaluations; 

• rj^: 2ksm flops; 

• M^: 2s^m flops; 

• A^: 4sm^ flops; 

• 7^+-"^: sm flops. 

Consequenly, this latter algorithm has a fixed computational cost of 1 Jacobian evaluation 
and flops, plus a cost of s function evaluations and 4sm^ + Aksm + sm flops per 
iteration. 
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Table 4.2: Outer-inner iteration for the blended implementation of HBVMs. 



n 
e 

for 



= I-ihQJo 

= /, (g) 



given 
i = 0,1,... 

f = fi/) 

v^ = -¥ + {vJ^)®if^ 

A^'O = 

for r = 0, 1, . . . 

if r > 



% actually, $7 is factored LU 
% e.g., f = 



[/, Jo] A"' 



Lr 



else 

w ' = r) 

end 



end =^ 
7^+1 = 7^ + 



returns 



end 



Table 4.3: Nonlinear iteration for the blended implementation of HBVMs. 





= I-ihQJo 




e 


= Is (g) 


1 % 


actu 




given 


% e.g 




for 


£ = 0,1,... 






= e<S 


) yo + {his) i 


5 /7 




f = fi/) 








If 










A' = e 


'e{u^ - rf) - 






f+i = 7^ + A^ 





end 



Chapter 5 

Line Integral Methods 



Sometimes conservative problems are not in Hamiltonian form and/or they posses multiple 
constants of motions, which are functionally independent. In certain cases, it is crucial to 
be able to preserve all of them, in order to obtain a faithfully simulation of the underlying 
dynamical system. For this reason, we extend the polynomial methods studied above, 
in order to cope with these, more general, conservative problems. The material of this 
chapter is based on [6l [5] . 



5.1 Introduction 

In the previous chapters, we have studied polynomial methods for approximately solving, 
on the interval [0,h], the initial value problem 



y'ich) = f{y{ch)) = J2^,{y)Pj{c), c G [0, 1], y{0) = yo e 
lM = I P,{r)fiyirh))dT, j>0. 



2m 



(5.1) 



(5.2) 



The methods that we have considered are characterized by a suitable polynomial a G n<j 
such that 

s-1 
j=0 



AC), 



ce[0,l], aiO)=yo, 



(5.3) 



then approximating the integrals {a) by a suitable quadrature formula. When the prob- 
lem is Hamiltonian, that is, /(•) = J'WH{-), with = —J = J~^, energy is conserved. 



for the discrete-time dynamical system defined by (5.3). Indeed, 



H{a{h)) - H{a{0)) = I V H{a{t)fa' {t)dt = h I VH{a{Th)fa'{Th)dT 

Jo Jo 

= h fvH{a{Th)fY.^,{a)P,{T)dT 



VH{a{Th))Pj{T)dT 



Uo 



s-1 

i=o 

s-1 
j=0 
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due to the fact that J is skew-symmetric. 

Let now consider the case where problem (5.1 ) is a general conservative (not necessarily 
Hamiltonian) problem, whose dimension will be denoted by m, for sake of brevity. Assume 
that it possess a set of u smooth, functionally independent (clearly, u < m), invariants. 
That is, there exists 

L: 

such that 

VL{yff{y) = 0, Wye R^, (5.4) 

where VL{y)'^ denotes the Jacobian matrix of L. In such a case, along the solution of 
(5.1 ), one has: 



dt 



L{y{t)) = WL{y{t)Yy'{t) = W L{y{t)Y f {y{t)) = Oe 



because of (5.4). Consequently, 

L(y(t))=L(yo), 



Vt > 0. 



From (5.1) and (5.5) one obtains: 



(5.5) 



L{y{h))-L{ym 



VL{y{t)ff{y{t))dt 
h l\L{y{ch)fY.P^^^)^^^y)^^ 

i>o 



where 7j(i/) is formally still defined by (5.2) and, moreover, we have set 



(5.6) 



P,{c)VL{y{ch))dc, J > 0. 



(5.7) 



Remark 5.1 It is important noticing that, because of (5.4), (5.6) continues to hold, if 
we replace y{t) by any other path (j{t). 



We observe that, because of the result of Lemma 2.2, one has (assuming, for sake of 
simplicity, that both L{y{t)) and f{y{t)) can be expanded in Taylor series at t = 0): 

0,(y) = O(/i^), j>0, (5.8) 



However, if we replace the infinite series in (5.1) with the finite sum in (5.3), one obtains, 
because of (5.4) and repeating similar steps as above: 

s-l 

L{a{h)) - L{am = hJ^M^f^A^) = -hY^M^f^A^) = 0{h'^^'), (5.9) 

j=0 j>s 



since (see (5.6)) 



5^0,(o-)%(a) = 0. 

j>0 



5.1. INTRODUCTION 
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In order to get conservation for the polynomial dynamical system, we perturb (5.3) as 
follows: 

s-l 

a\ch) = Y,l,i^)Pjic) - M^)(^, ce[0,l], a(0) = yo, (5.10) 

j=0 



where 0o('^) is defined according to (5.7), and a G M'^ is determined in order to enforce 



the conservation of the invariants. By repeating similar steps as above, we obtain: 

= L{a{h)) - L{a{0)) = h[ VL{a{ch)fa'{ch)dc (5.11) 

Jo 

s-l 

= h^(j)j{af-fj{a) - h[(t)o{af(t)o{(x)]a. 

j=0 

Consequently, conservation is gained, provided that 



s-l 



o(a)^0o(a)]« = 5^0,(a)% (5.12) 

j=0 

The following result holds true. 

Theorem 5.1 The vector a exists and is unique, for all sufficiently small step sizes h 
and, moreover, 

a = 0{h^'). (5.13) 

Proof. If the invariants are functionally independent, cr(0) is a regular point for the con- 
straints, so that VL((t(0)) has full column rank (i.e., u). Considering that 

0o(a) = [ VL{a{ch))dc -> VL{yo), as /i -> 0, 
Jo 

one has that matrix 

Mo = [M^^fM^)] 

is symmetric and positive definite and, therefore, nonsingular. The existence and unique- 



ness of a then follows from the Implicit Function Theorem. Moreover, since (see (5.8)) 

Mo = 0(/i°), 

then 

s-l 

a = Mo"^5^0,(a)%(a) = -Mo'^ 0,(a)^,(a) = 0{h'^). 

j=0 j>s 

This completes the proof. □ 



We now consider the following question: i.e., the polynomial a as defined in (5.3) 



doesn't satisfy, in general, L{a{h)) = L^yo), even though a{h) — y{h) = 0{h'^'^^^). Con- 



versely, for the polynomial a defined by (5.10)-(5.12) one has L{a{h)) = L{yo). Moreover 



next theorem states that its order of convergence remains the same. 

Theorem 5.2 Let a be defined according to (f5.iO|j-(f5J^. Then a{h)—y{h) = 0(/i^''+^). 
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Proof. By using similar steps as those used in the proof of Theorem 3.1, one has: 



a{h)-y{h) = y{h;h,a{h))-y{h;0,a{0)) 



dt 



y{h;t,a{t))dt 



(^y{h;9,a{t)) + ^y(h;t,co) (r'{t)]dt 

\00 e=t OOO uj=a{t) 

%(/i,t)[-/(a(t)) + a'(t)]dt 





hf ^{h,Th)[-f{o{Th))+a'{Th)]dT 

Jo 



-h / <l>{h,Th) 



s-l 



Li>o 



j=Q 



dr 



-h [ <l>{h,Th)^Pj{T)jj{a)dr ~ h [ <l>(/i, r/i)0o(fx)adr 



= G{Th) 



^{h,Th) Pj(r)dr 



— V 

--0(,hi) 



= 0(1) 



j>s 



5.2 Discretization 



As was previously observed, (5.10)-(5.12) doesn't yet define a method, but a conservative 



formula. As matter of fact, a numerical method is obtained when the integrals (see (5.2) 
and ( [K7| )) 

7j(cr), 0j(cr), j = 0,...,s-l, 

are approximated by means of a suitable quadrature formula. In principle, they could be 
approximated by means of different quadrature formulae: 

• one formula, based at the abscissae < ci < . . . < < 1 and corresponding weights 
{bi}, of order q, for approximating 7j(o"): 



i = 0, . . . , s — 1, 



with 



A,(/i) = 0(/i^-^), j = 0,...,s-l; 



(5.14) 



(5.15) 



another formula, based at the abscissae < ri < . . . < < 1 and corresponding 
weights {A}; of order q, for approximating (pj{cr): 



i = 0, . . . , s — 1, 



with 



j=0,...,s-l. 



(5.16) 
(5.17) 
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In the following, for sake of simplicity, we shall consider the following choices of such 
abscissae: 

Pk{c^) = 0, t = l,...,k, q = 2k, (5.18) 

Prir,) = 0, j = l,...,r, ^ q = 2r. (5.19) 



Definition 5.1 We shall refer to such a method as LIM(r, k, s), where LIM is the acronym 
for Line Integral Method. 

Remark 5.2 We observe that: 

• LIM{0, s, s) is the s-stage Gauss method, 

• LIM{0,k,s) is the HBVM{k,s) method, 

where r = means that no invariant conservation is seeked. 

After discretization, the polynomial a is obviously formally replaced by the polynomial 
M G n, such that: 



s-1 



u 



\ch) = Y,l,P'Ac) - <^o«, CG [0,1], 



m(0) = ?/o, 



with, in general, (see (Q, ^Jj, ( |5.14D -( |5.17D ) 

k 

Ij = ^biPjici)f{u{cih)) = -fj{u) + Aj{h), 

i=l 
r 



(5.20) 



(5.21) 
(5.22) 



for j = 0, . . . , s - 1, and (see (|518|)-(|539 

s-1 



s-1 

^ = ll^hj = 5^[0,(n) + vl/,(/i)r[7,(n) + A,(/.)] 

j=0 j=0 



s-1 

E 

j=0 



--0(h^^) 



= 0(/i2fc) 



--0{h' 



2{r+fc-j) 



0,(«)^ 7,(«).+ ,^,(/^)'^, 7.(«) + 0.W +^,ihfA,ih) 

= 0(h2fc-i) 



= 0{/iJ) =0{hJ) =0(/i2'-i) 
s-1 

j>s j=0 

0{h^') + Oih^') + 0{h^'') + 0(/i2('^+'=-^+i)). (5.23) 



Remark 5.3 From (5.21)-(5.22), and (5.23), it follows that the (block) size of the discrete 
problem is 2s + 1.- 

• the s coefficients 7j G M™", j = 0, . . . , s — 1, 
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• the s coefficients (f)j G W^^" , j = 0, . . . , s — 1, 

• the vector d G M'^. 

even though it must be stressed that (pj is a matrix with v columns, whereas jj is a vector. 
The efficient implementation of such methods is, however, still under investigation. 

The following result then holds true. 

Theorem 5.3 For all r > s and k > s, for LIM{r, k, s) one obtains: 

u{h)-y{h) = 0{h^'+'). 



Proof. From (5.23) and the hypotheses r > s and k > s, it follows that 

a = 0{h^'). 



Moreover, ones has, by repeating similar steps as in Theorem 5.2 



u{h)-y{h) = y{h-h,u{h))-y{h-Q,u{Q)) 



dt 



y{h-Xu{t))dt 



^ y{h-e,u{t)) + -^y(h;t,tu) uit)]dt 



09' 



dco' 



uj=u{t) 



^h,t)[-f{u{t))+u\t)]dt 



h $(/i,r/i)[-/(n(r/i)) + n'(r/i)]dr 



-h / <l>{h,Th) 



.i>o j=o 



dr 



h / <i>{h,Th) 



.j>o 



s-1 



3=0 



= 0(/i2fc-i)' 

7i(w)- 



= 0(h2'-)' 

0o(m) - ^o(/i) 



-h / ^{h,Th)y2Pj{r)jj{u)dT - h / ^{h,Th)(f)o{u)adT 
Jo j.>, Jo 

pi s-l pi 

-h / $(/i,r/i) VP,(r)Aj(/i)dr + /i / ^{h,Th)^o{u)adT 

Jo ■_r, Jo 



a 



dr 



s-1 
j=0 



^{h,Th)Pj{T)dT 

V ' 

= 0{hi) 



(!>{h,Th)Pj{T)dT 



— V 

--0{h3) 



= 0{hi) 

7j(m) -h 



= 0(/l2'=-J) 



$(/i,r/i)dr 



= 0(1) 



^{h, Th)dT 



'oil)' 



= 0(ft2») 

= o(i) 

= 0(/l2r) 
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Consequently, the following statement easily follows. 

Corollary 5.1 For all r > s and k > s, LIM{r, k, s) has order 2s. 

Concerning the conservations of the invariants, the following result holds true. 

Theorem 5.4 For given r,k > s, LIM{r, k, s) exactly conserves polynomial invariants of 
degree 

(5.24) 
(5.25) 



z/ < 



2r 



For all suitably regular invariants, 

L{u{h))-L{yo) = 0{h''^+^). 



Proof. One has, by virtue of (5.20)-(5.22): 



L{u{h)) - L{yo) = L{u{h)) - L{a{0)) = h I VL{u{ch)fu{ch)dc 



s-1 



j=0 

= El. 



In case L is a polynomial of degree u satisfying (5.24), the quadrature formula (5.22) is 
exact, so that 



J = 0,...,s-1. 



Consequently, El = 0, by virtue of (5.23). This proves the first part of the thesis. In any 
other case, one has: 



Er 



h 



s-1 



^(0j-^j)% - ((0o-^oJ^</>o) a 

.j=0 



s-1 s-1 

J]0j7, - (00 0o) a - ^Jlj + (£o^ 



a 



j=0 



i=o 



=0, see ( |5.23| ) 
thus proving ( 5.25[ ). □ 



Remark 5.4 From (5.25), one obtains that, for any suitably regular set of invariants. 



conservation can always be practically obtained, provided that r is large enough. Indeed, 
it is enough to obtain conservation up to roundoff errors. 

Moreover, if some of the invariants are polynomials of low degree, then, in principle, 
a less accurate quadrature formula could be used for approximating the corresponding 
integrals (5 
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The following result can be also proved, by using arguments similar to those used in 
Section 13.41 



Theorem 5.5 Provided that the abscissae (5.18)-(5.19) are symmetrically distributed in 



the interval [0,1], LIM{r\k^ s) is a symmetric method. 

We now provide a couple of straightforward fully-conserving generalizations of HBVMs 
and Gauss-Legendre Runge-Kutta methods. 

5.2.1 LIM(A;, A;,s) 

Such conserving methods can be regarded as a straightforward generalization of HBVM(/c, s) 
methods. In such the same set of abscissae, 

< ci < . . . < Cfc < 1, Pfc(ci) = 0, z = l,...,/c, 

are used for the quadraures approximating 7j(tt), (pj{u), j = 0, . . . , s — 1. Consequently, 
by setting as usual Yi = u{cih), one obtains: 



s-l 



k 



i — 1, . . . , fc. 



J = 0, . . . ,s - 1, 

k 

Y,biP,{c,)VL{u{Ye)), 



£=1 



0^00 1 a 



j=0 



with the new approximation given by 



(5.26) 



yi=yo + h^ hif{Yi) - /i0oa- 



i=l 



5.2.2 LIM(fc,s,s) 



(5.27) 



These methods turn out to be fully conservative variants of the s-stage Gauss methods. 
In such two sets of abscissae are used, i.e.. 



< Cl < . . . < Cs < 1, Ps{ci) = 0, i = 1, . . . , s 

< Ti < ... < Tfc < 1, Pfc(r,) = l,...,s. 
The resulting method can be easily seen to be, by setting Yi = u{cih): 

Yi = Vo + hy^ Pjix)dx% - h^^a, 

s 

Zi = '^Lis{Te)Yi, 

i=l 



i = l,...,s, 



1, . . . , fc. 



(5.28) 
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% = 5^6,p,(Q)/(r,), 

£=1 

j = 0,...,s-l, (5.29) 

k 

0, = ^/3,P,-(q)VL(Z, 



e=i 

i=o 



where the Lis{T) are the Lagrange polynomials defined at the abscissae (5.28), and the 
new approximation is given by 



s 



yi = yo + hJ2 ^^10"^) - hka. (5.30) 

i=l 

In this case, the first equation in ( |5.29[ ) can be recast in the more usual form, 

u{cih) = f{u{c.ih)) - 0o«, i = 1, . . . , fc, 

which emphasizes the connection with the collocation conditions of the original s-stage 
Gauss method. 

5.3 Numerical tests 

In this section, we report a few numerical tests on a couple of conservative problems, 
possessing multiple invariants. 

The Kepler problem 

A noticeable example of Hamiltonian problem with multiple invariants, is the Kepler 
problem, defined by the (non-polynomial) Hamiltonian: 

^(q,P) = J||P||2-^, q,pGM2, (5.31) 
^ liq||2 

It admits the following invariants of motions, besides the Hamiltonian: 

• the angular momentum, 

^(q, P) = (I1P2 - q2Pi (5.32) 

which is a quadratic invariant; 

• the so called Laplace- Runge-Lenz (LRL) vector which, for the problem at hand, 
implies the conservation of the following quantity, 

^(q, p) = <i2vl - qivm - ttV- (5.33) 

llq||2 

When starting at the initial point 



gi = 1 - £, g2 = Pi = 0, P2 = Y ^37' (^-^^^ 
it has a periodic orbit of period T = 27r, which is, in the (gi, g2)-plane, an ellipse of 



eccentricity e. We first consider the integration of problem (5.34) with e = 0.6, by using 
the following fourth order methods with a constant stepsize h = lO^^vr: 
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• LIM(0,2,2), (i.e., the 2-stage Gauss method); 

• LIM(0,8,2), (i.e., HBVM(8,2)); 

• LIM(8,2,2), (i.e., the "fully conservative" variant of the 2-stage Gauss method); 

• LIM(8,8,2), (i.,e., the "fully conservative" variant of the HBVM(8,2) method). 

Indeed, for the "fully conserving methods" , r = 8 is enough to obtain a practical conser- 
vation of all invariants. 



In Figures 5.1 and 5.2 is the plot of the error in the invariants for LIM(0,2,2) ans 
LIM(0,8,2), respectively: for the first method, the error in the Hamiltonian is bounded, 
the angular momentum (which is a quadratic invariant) is conserved up to roundoff, and 
a drift apparently occurs in the LRL invariant; for the second method the situation is 
similar, with the roles of the Hamiltonian and of the angular momentum exchanged each 
other. Also in this drift seems to occur in the LRL invariant: this is better 



evidenced in Figure 5.3, where the drift of the LRL invariant (5.33) is plotted over a 
longer interval, and it appears to be same for both LIM(0,2,2) and LIM(0,8,2). 

On the other hand, both LIM(8,2,2) and LIM(8,8,2) conserve all the invariants up to 



roundoff. In Figure 5.4, there is the plot of the error in the numerical solution (measured 
at each period) for all methods: it is evident that its growth is linear, and with the 
same order, even though LIM(0,2,2) (i.e., the two stage, symplectic Gauss method) is less 
accurate than the other methods. 

This scenario changes when the eccentricity e ^ 1: indeed, in such constant 
stepsize is very inefficient and a variable stepsize would be preferable. By using a standard 
mesh selection strategy, based on the control of the local error, e.g., 

1 

(tol \ 

where: 

• hold is the current stepsize; 

• knew is the new stepsize; 

• 0.85 is a "safety" factor; 

• tol is the prescribed tolerance for the local error; 

• e is an esatimate of the latter error; 

• p is the order of the method; 

it is known that symplectic methods may suffer from a quadratic error growth (instead 
of a linear one, as shown, e.g., in [iQl pp. 303-305]). Let us then see what happens when 
e = 0.99, and a tolerance tol = 10~^ is used for all the above methods: 



for the 2-stage Gauss method, from the plot in Figure [5~5] one has that a drift in both 
the Hamiltonian and the LRL invariant is present, whereas the angular momentum 
is preserved up to roundoff; 



for the HBVM(8,2) method, from the plot in Figure 5.6 one has that a drift in both 
the angular and the LRL invariant is present, whereas the Hamiltonian is preserved 
up to roundoff; 
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t 



Figure 5.1: Kepler problem, e = 0.6, 2-stage Gauss method, h = W ^vr, invariants errors. 




Hamiltonian 
Angular momentum 



100 200 300 400 500 600 

t 

Figure 5.2: Kepler problem, e = 0.6, HBVM(8,2) method, h = lO^^vr, invariants errors. 
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t 



Figure 5.3: Kepler problem, e = 0.6, drift in the LRL invariant for both the 2-stage Gauss 
method and the HBVM(8,2) method, h = lO^^vr. 




Figure 5.4: Kepler problem, e = 0.6, error in the numerical solution, h = 10 ^vr. 
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for both LIM(8,2,2) and LIM(8,8,2), all the invariants are conserved up to roundoff. 



At last, in Figure [57f| there is the plot of the error, measured at each period, for 100 peri- 
ods: it is evident that the solution provided by the 2-stage Gauss method soon becomes 
meaningless, whereas for all the other methods a linear error growth is observed, the fully 
conserving methods being slightly more accurate than HBVM(8,2). 

The Lotka-Volterra problem 

The second test problem that we consider is the Lotka-Volterra problem, i.e., a problem 
in Poisson form, 

y' = B{y)VH{y), y(0) = yo, (5.35) 

with 

B{yf = -Biy), Vy, 

and the scalar function H{y) is still called the Hamiltonian. Also in this case, the Hamil- 
tonian is a constant of motion, since: 

^^H{y{t)) = VH{y{t)fy'{t) = VH{yfB{ymH{y{t)) = 0, 

B{y(t)) being skew-symmetric. Moreover, each function C(y) such that 

\/C{yfB{y) = 0, (5.36) 
is an invariant for the corresponding dynamical system. Indeed, one has: 

^C{y{t)) = VC{y{t)fy'{t) = VC{yfB{y{t))VH{y{t)) = 0, 



because of (5.36). A function C{y) which satisfies (5.36) is said a Casimir ioi (5.35). Let 
us then consider the following problem [32j, for which y = (|/i, ?/2, l/s)^, 

/ cyiy2 bcym \ 
^(y) = -cyiy2 -7/21/3 , abc = -l, 
\ -bcym ?/2l/3 J 

the Hamiltonian is 

H{y) = abyi + y2 - ay^ + log ?/2 - log 2/3, 
and, moreover, there is the following Casimir: 

C(y) = a6 log 1/1 - 61og?/2 + logi/3. 
By using the following set of parameters, 

a = —2, b = —1, c = —0.5, = 1, = 2, 

and initial point: 

yo = ( 1 1.9 0.5 
the solution turns out to be periodic with period 

T ^ 2.878130103817. 
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Figure 5.5: Kepler problem, e = 0.99, 2-stage Gauss method, tol = 10 ^, invariants errors. 
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Figure 5.6: Kepler problem, e = 0.99, HBVM(8,2) method, tol = 10 ^, invariants errors. 
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Figure 5.7: Kepler problem, e = 0.99, error in the numerical solution, tol = 10 
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We now solve this problem with a constant stepsize 

h = T/30 ^ 0.096, 

so that we can check both the errors in the solution and in the invariants. We solve, at 



first, the problem by using the 2-stage Gauss method. In Figure 5.8 we plot the errors in 
the invariants along the numerical solution: as one can see, both of them exhibit a drift. 
Then, we use LIM(8,2,2), with the same stepsize, by imposing only the Hamiltonian 
conservation: indeed, in the case of the Kepler problem, this was sufficient to obtain a 
linear growth for the solution error. In Figure 5.9| we plot the error in the numerical 



Hamiltonian and Casimir, thus showing a practical conservation of the former, and a 
linear drift for the latter. Finally, we use LIM(8,2,2), with the same stepsize, by imposing 
both the conservation of the Hamiltonian and of the Casimir, which are conserved up to 



roundoff. At last, in Figure 5.10 we plot the measured error in the solution, measured 
over 100 periods: one then concludes that a linear error growth is attained only when 
preserving both the invariants; conversely, a quadratic error growth occurs. 
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Figure 5.8: Lotka-Volterra problem, 2-stage Gauss method, h = T/30, invariants errors. 
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;;ure 5.9: Lotka-Volterra problem, LIM(8,2,2) method with only the Hamiltonian preserved, 
= T/30, invariants errors. 
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Figure 5.10: Lotka-Volterra problem, solution errors with stepsize h = T/30. 



Chapter 6 

Further developments and references 



We end these lecture notes, by adding that further interesting developments, such as the 
possibility of getting, in a weakened sense, methods which are both symplectic and energy- 
conserving, have been considered in |20l[T5]. We also mention that a noticeable extension 
of this approach, for PRK methods, has been recently devised in [60]. A further line of 
investigation deals with multistep energy-preserving method, as is sketched in |19]. Last, 
but not least, the efficient implementation of such methods deserves to be investigated as 
well. 
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